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The numerical study of anyonic systems is known to be highly challenging due to their non- 
bosonic, non-fermionic particle exchange statistics, and with the exception of certain models for 
which analytical solutions exist, very little is known about their collective behaviour as a result. 
Meanwhile, the density matrix renormalisation group (DMRG) algorithm is an exceptionally pow¬ 
erful numerical technique for calculating the ground state of a low-dimensional lattice Hamiltonian, 
and has been applied to the study of bosonic, fermionic, and group-symmetric systems. The re¬ 
cent development of a tensor network formulation for anyonic systems opened up the possibility of 
studying these systems using algorithms such as DMRG, though this has proved challenging both in 
terms of programming complexity and computational cost. This paper presents the implementation 
of DMRG for finite anyonic systems, including a detailed scheme for the implementation of anyonic 
tensors with optimal scaling of computational cost. The anyonic DMRG algorithm is demonstrated 
by calculating the ground state energy of the Golden Chain, which has become the benchmark sys¬ 
tem for the numerical study of anyons, and is shown to produce results comparable to those of the 
anyonic TEBD algorithm and superior to the variationally optimised anyonic MERA, at far lesser 
computational cost. 


I. INTRODUCTION 

The Density Matrix Renormalisation Group (DMRG) 
algorithm 1-3 is an algorithm for computing the ground 
state of a low-dimensional lattice Hamiltonian, and is 
arguably one of the most successful and widespread al¬ 
gorithms in condensed matter physics. 3 With the advent 
of tensor network representations of quantum states, it 
is now recognised that DMRG may be understood as a 
variational tensor network algorithm to compute a Ma¬ 
trix Product State (MPS) representation of the ground 
state of a quantum system. 4,5 

Anyons, meanwhile, are quasiparticles capable of ex¬ 
isting only in one- or two-dimensional systems, 6 and 
are extremely challenging to simulate numerically due to 
their non-bosonic, non-fermionic exchange statistics. To 
date, simulation of anyonic systems has largely consisted 
of exact diagonalisation performed on chains or ladders 
of around 30-40 lattice sites,' -12 and exact calculations 
where mappings to analytically solvable models exist. 7,13 
Where available, such mappings have also been exploited 
to permit some systems with quantum group statistics 
to be studied using DMRG. 14-16 The study of anyonic 
systems has so far been largely theoretical, though the 
recently reported detection of Ising anyons at the ends of 
iron nanowires 1 ' may be about to change all that, lending 
pressing urgency to the development of efficient numeri¬ 
cal algorithms for the simulation of anyonic systems. It 
would therefore be highly advantageous to convert exist¬ 
ing, highly effective algorithms such as DMRG for use 
with arbitrary anyonic systems. 

With this goal in mind, recent work on symmetries 
in tensor network states 18 led to the development of a 
tensor network formalism for anyonic systems. 19,20 This 


formalism permits existing tensor network algorithms to 
be adapted to the study of anyonic systems, including 
the scale-invariant Multi-Scale Entanglement Renormal¬ 
isation Ansatz (MERA) 19,21-23 and the Time-Evolving 
Block Decimation (TEBD) algorithm. 24-26 In the present 
paper, we show—in quite considerable detail how the 
DMRG algorithm may also be implemented for anyonic 
systems. Our aims in doing so are twofold: First, how 
could we resist bringing together the most challenging 
particles to simulate, and the most powerful numerical 
technique for low-dimensional systems, especially when 
those particles can only exist in one or two dimensions? 
Second, the DMRG algorithm is the perfect context in 
which to provide a more detailed exposition of the any¬ 
onic tensor formalism first introduced in Ref. 19 and its 
application. Its tensor networks are more complex than 
those of TEBD but less so than those of the MERA, and 
so it is possible to present the algorithm at a level of de¬ 
tail which is simply not practical for the MERA, while 
nevertheless involving a level of sophistication which is 
not required for the simpler tensor networks encountered 
in TEBD. 

This paper is organised as follows: 

In Sec. IIA we discuss the construction of anyonic 
tensors, the choice of bases on tensor indices, and how 
these relate to fusion tree diagrams. In Sec. IIB we see 
how these tensors may be efficiently manipulated, and in 
Sec. IIC we construct anyonic tensor networks, examine 
the pairwise contraction of anyonic tensors, and compare 
it with the whole-network approach employed in Ref. 26. 
We also compare and contrast the anyonic tensor net¬ 
work formalism with that used for conventional tensor 
networks. In Sec. Ill A we use anyonic tensors to de¬ 
fine an MPS Ansatz for anyonic systems on the disc, and 
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in Sec. Ill we show how the anyonic version of the fi¬ 
nite DMRG algorithm uses this Ansatz to compute the 
ground state of a Hamiltonian on a finite chain on the 
disc. Section IV briefly addresses the extension of this 
Ansatz to the study of chains on the torus as an example 
of anyonic DMRG on surfaces of higher genus. 

The infinite DMRG algorithm is not addressed in this 
paper: Its adaptation to anyonic systems is less straight¬ 
forward than that of the finite DMRG algorithm, and 
caution is required to avoid inadvertently imposing non¬ 
physical constraints on the states which may be explored. 
The nature of these constraints and techniques for their 
avoidance are discussed in Ref. 27, with their application 
to anyonic systems being addressed explicitly in Sec. Ill F 
of that paper. 


II. ANYONIC TENSORS 

This Section details the construction of anyonic ten¬ 
sors and tensor networks. A normalisation scheme is 
introduced which is compatible with the diagrammatic 
isotopy convention described in Refs. 28 and 29. Famil¬ 
iarity with the diagrammatic representation of anyonic 
states employed in these works is assumed. The focus 
of this Section is pragmatic, with emphasis on concerns 
such as choices of normalisation convention, and is in¬ 
tended to complement the more rigorous but also more 
abstract treatment presented in Ref. 19. 


A. Construction 

1. Anyonic state vectors 

As in Ref. 19 we begin by considering a surface of 
genus 0 on which there exists a lattice C of n sites, pop¬ 
ulated by anyons [Fig. 1 (i)]. We then choose a basis by 
introducing a linear ordering of these n sites, as discussed 
in Ref. 30 and illustrated by the dashed line in Fig. l(i), 
and choosing a unidirectional 31 fusion tree [Fig. l(ii)], 
whose projection onto the manifold corresponds to the 
linear ordering adopted in Fig. l(i). Extension of the fu¬ 
sion tree to the boundary of the disc permits inclusion of 
a boundary charge a&. This charge may be represented 
in two ways: As a total charge at the bottom of the fu¬ 
sion tree, or as a leaf carrying charge ajj. We favour the 
latter approach, as a fusion tree may then be understood 
as a recipe for the construction of a state, the lines of 
the tree corresponding to world lines, with cross-sections 
through the tree being edge-on views of the manifold at 
a given instant in time. 32 Given the equivalence between 
a boundary charge and an additional lattice site, we will 
not need to explicitly consider the possibility of a bound¬ 
ary charge beyond this point. 

A state is defined as a weighted sum over labellings 
of the fusion tree which are consistent with the anyonic 



FIG. 1. Choosing a basis for a system of four anyons on 
the disc (plus a possible boundary charge): (i) Four anyons 
on the disc. A linear ordering, represented by the dashed 
line, has been imposed, (ii) A fusion tree is selected whose 
projection onto the disc coincides with the imposed linear 
ordering. The labels on this fusion tree represent a selection 
of combined charges. For example, ay is the total charge of 
sites 1 and 2 together, and ag is the total charge of sites 1, 2, 
and 3 together, (iii) A state is constructed as a weighted sum 
over valid labelings of the fusion tree. 


fusion rules, and a total charge of I [Fig. 1 (iii)] - The num¬ 
ber of indices involved may be reduced considerably by 
introducing an index fij which enumerates all labellings 
of the fusion tree consistent with these conditions. Any 
labelling of a fusion tree may then be identified with a 
pair (a, fi a ) where a is the total charge of the fusion tree, 
and the admissible values of index [i a enumerate the la¬ 
bellings of the tree consistent with total charge a. (For 
example, if there are seven valid labellings with a given 
value of total charge a, then /r Q takes an integer value in 
the range [1,7].) When describing a state we will only 
need pairs for which a = I, but other values of a will 
be necessary when describing operators and three-index 
tensors. 

Using this more compact notation the coefficients de¬ 
scribing a state may now be written 

c (a,Ma), 

and if we introduce a linear ordering of the pairs (a,/r a ) 
which is enumerated by greek index a, then we obtain the 
state vector c a where a enumerates the basis elements 
corresponding to all valid labellings of the chosen fusion 
tree. 

For tensor network algorithms, it is conventional to 
represent a single-index tensor such as c“ by a circle (de¬ 
noting the tensor) with a single line emerging (denoting 
the index), as shown in Fig. 2(i). Ignoring, for the mo¬ 
ment, that a corresponds to labellings of a fusion tree, 
and writing the basis of our Hilbert space as |ii),..., |i a ), 
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FIG. 2. (i) Conventional tensor network notation for a vec¬ 
tor c“ (i.e. a single-index tensor), or associated state \if>) = 
5^Q Ca |*“)- The circle represents the tensor c, and the sin¬ 
gle emerging line represents the index. The conjugate state 
(-01, and associated tensor 4, is represented by diagram (ii). 
(iii) Diagram in conventional tensor network notation corre¬ 
sponding to the inner product c“cj, (repeated index summed). 
The line connecting c and c4 represents a summed index ap¬ 
pearing on both tensors, (iv) Extension of graphical notation 
for an anyonic vector c“, incorporating the fusion tree whose 
labellings are enumerated by a. The fusion tree is drawn in 
grey to indicate a normalisation convention where no numer¬ 
ical factors are associated with vertices or loops, (v) The 
graphical representation of the conjugate state is obtained by 
vertical reflection and reversal of arrows, (vi) The anyonic in¬ 
ner product, which reduces in this normalisation convention 
to diagram (vii), again corresponding to c a c' a . 


we may identify a vector c“ with a state 

\^) = Y J C a \i a ). ( 2 ) 

a. 

Its conjugate, 

M = £4<i“l, ( 3 ) 

a 

is obtained by taking the Hermitian conjugate of c a to 
obtain c' a , i.e. 

4 = (c“)* V a, (4) 

vertically reflecting Fig. 2(i), and reversing all arrows 
(which is equivalent to conjugation of charges), to ob¬ 
tain Fig. 2(ii). Note that we write indices upper or lower 
to match the upgoing or downgoing orientations of the 
associated legs in the graphical notation. It is convenient 
to have the vector representation of a physical state sat¬ 
isfy the normalisation condition 

W#) = c“4 = 1 (5) 

(where there is an implicit sum over the repeated index 
a), and this is represented diagrammatically as shown in 
Fig. 2 (iii). 

For anyons, as shown in Fig. 2(iv), we add in a dia¬ 
grammatic representation of the fusion tree by drawing 




FIG. 3. Example diagrammatic representations of (i) an any¬ 
onic operator and (ii) an anyonic three-index tensor (this 
example has two downgoing indices and one upgoing index; 
one downgoing index and two upgoing indices is also accept¬ 
able). Note that the tensors M and T implicitly play the 
role of vertices connecting the upper and lower fusion trees, 
and thus matrix Mg contains non-zero entries only for values 
of a = {a, Ha) and /I = (b, fib) such that a = b. Similarly, 
if 7 = (c, He) then T'“ /3 is non-zero only where the product 
a x b —> c is permitted by the fusion rules of the anyon model. 


it on the open end of the index, illustrating the basis 
whose elements are enumerated by a. In this paper we 
draw the fusion tree in pale grey to indicate that this tree 
is not associated with any diagrammatic factors as per 
Refs. 28 and 29. The conjugate of Fig. 2(iv) is shown in 
Fig. 2(v), and is again obtained by vertical reflection of 
the diagram and complex conjugation of the coefficients 
c“ to obtain c' a . Finally, Fig. 2(vi) is the anyonic ver¬ 
sion of the inner product, Fig. 2(iii). To evaluate this 
diagram, first confirm that the two pale grey fusion trees 
are mirror images of one another, indicating that c“ and 
4 are in conjugate bases, and then eliminate them to ob¬ 
tain Fig. 2(vii) [which is simply Fig. 2(iii) supplemented 
by an orientation arrow on the shared index] and thus 
Eq. (5). Note that no numerical factor arises from the 
elimination of the grey fusion trees, as they merely indi¬ 
cate the basis employed on index a. Because these fusion 
trees are not associated with the diagrammatic equations 
and normalisation factors of the convention described in 
Refs. 28 and 29, we call these implicit fusion trees, being 
implied by the choice of basis on a. We will contrast 
this with explicit fusion trees, which are associated with 
the numerical coefficients and rules of the diagrammatic 
isotopy convention, in Sec. IIB 3. 


2. Two-index operators and three-index tensors 

This notation generalises to operators (matrices) and 
three-index tensors as shown in Fig. 3. Note that for ob¬ 
jects with two or more indices, the total anyonic charge 
on a fusion tree may take values other than I. That is, 
if an index a is decomposed into a pair (a, Ha), then all 
values of a are allowed provided they are associated with 
valid labellings of the fusion tree. This is particularly im¬ 
portant for operators acting on a subsystem, which may 
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have non-trivial total charge even though the system as 
a whole carries a total charge of I. We find it convenient 
to adopt a convention where all arrows on fusion trees 
are upgoing, and the entries of a matrix Mg = M^’^) 
are then constrained to be non-zero only for a = b. For a 
three-index tensor an entry may be non-zero only where 
the implicit vertex between a , b, and c is permitted by 
the fusion rules. 


3. Charge multiplicities 

Considering a three-index tensor TF as in Fig. 3, sup¬ 
pose that the fusion rules admit outcomes with multiplic¬ 
ity greater than one, for example 


These include hermitian conjugation of multi-index ob¬ 
jects, fusing and splitting of indices, conversion between 
upper and lower indices using caps and cups, and changes 
of fusion tree basis using braiding and F-moves. 


1. Hermitian conjugation 

The Hermitian conjugation of anyonic vectors in 
the graphical notation was introduced in Eq. (3) and 
Fig. 2(v). Extension to multiple-index objects (for ex¬ 
ample T“^) is straightforward: All entries in the object 
are replaced by their complex conjugates, while all upper 
indices become lower indices and vice versa while retain¬ 
ing their left-to-right ordering, e.g. 


8 x 8 ->l + 8 + 8 +.... ( 6 ) 

The identical fusion outcomes are commonly distin¬ 
guished by means of a multiplicity index associated with 
the fusion vertex, denoted /r. In the graphical represen¬ 
tation, however, a three-charge vertex carries only three 
indices. Consequently, one of the charges iterates not 
only over the output charges, but also over the associ¬ 
ated degeneracy index. Thus, given the fusion rule in 
Eq. ( 6 ) and component a of a and component b of fj 
satisfy a = b = 8 , then component c of label 7 takes 
values in { 1 , 81 , 82 ,...} with the subscript denoting the 
multiplicity index. 

If values in a and b also carry multiplicity indices, then 
c may carry as many different multiplicity indices as are 
required to distinguish each unique pair ab and fusion 
outcome. For example, if a G { 81 , 82 } and b G { 81 , 82 } 
then c G {li, I 2 , 13 , 14 , 81 , ..., 88 ,. . .}. There may come 
a point, however, at which information about the degen¬ 
eracy indices associated with the formation of charges a 
and b is no longer required. At this time, charges a and 
b may be collapsed such that 
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a G {8}, Ha= Hat 

2=1 

(7) 

2 

b G {8}, Hb = ^2Fbi, 

2=1 

(8) 


and similarly, 

CG { 81 , 82 }, n Cl =^2n Ci , fi C2 = ^ Ha (9) 

i odd i even 

where it has been assumed that odd i in Ha are associated 
with the first 8 in Eq. ( 6 ) and even i are associated with 
the second. 


B. Manipulation 

In this Section we present all single-tensor opera¬ 
tions required for the implementation of anyonic DMRG. 


^)l 0 = (Tf)* V a, /?, 7 , (10) 

and the associated graphical representation is reflected 
vertically. 


2. F-moves 


The simplest manipulation of a single tensor with 
which we will be concerned is the F-move, which is a 
change of basis on a composite index a corresponding to 
an alternative choice of implicit fusion tree. The F-move 
is a unitary transformation which may be applied to any 
part of a unidirectional fusion tree (one made entirely out 
of the same orientation of vertex) according to the rule 



where the ten-index tensor [Fd bc ](eap)(f is specified by 
the anyon model. In this paper we use only F-moves per¬ 
formed at the trunk of the fusion tree, though they may 
be freely performed at any level of the tree, on implicit 
or explicit indices. Note that for more extensive fusion 
trees, such as that shown in Figures l(iii) and 2(iv), it 
is common for multiple values of a to transform iden¬ 
tically. In this example, if performing a transformation 
at the trunk of the fusion tree, the F-coefficients are in¬ 
sensitive to the labels < 21 , a 2, <23, and < 27 , and it may 
be advantageous to treat these indices collectively as if 
they corresponded to a degeneracy of charge as- This 
is particularly true if the leaves of the fusion tree are 
supplemented with auxiliary degrees of freedom such as 
position data, or accessory qubits, as these may always 
be treated as degeneracies of the corresponding charge 
label, and thus all acquire the same multiplicative coef¬ 
ficients from operations on the fusion tree, consistently 
across multiple anyonic manipulations. 
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For inverse F-moves the relevant coefficients may be 
obtained by recognising that for fixed a, b, c, d the F ten¬ 
sor may be written as a matrix mapping between the 
initial and final bases, and this matrix is invertible, 

r Isabel 

L d J (e, a,P)(f,ti,v) L d j (f,!i,v)(g,p,a) 

(f,w) ( 12 ) 

= fi(e,a,(3)(g,p,cr) = ^eg^ap^(3a • 

For F-moves on the trees associated with kets, the coef¬ 
ficients may be obtained by vertical reflection, arrow re¬ 
versal (which takes the place of conjugation of charges), 
and complex conjugation of the coefficients of F or F -1 
as appropriate. 

3. Combining and splitting indices 

a. Index fusion and splitting operations: Consider 
again a vector c“ describing the state of five anyons on 
the disc, as per Fig. 2(iv). Suppose there exists an oper¬ 
ator Mfj acting on only three consecutive anyons, which 
we wish to apply to ai, 02 , and a 3 . We proceed as follows: 
First, we use F-moves to place c“ into a more convenient 
basis as shown in Fig. 4(i). We then introduce an explicit 
vertex normalised in accordance with the diagrammatic 
isotopy convention [Fig. 4(ii)]. Unlike the implicit ver¬ 
tices used so far, which merely provide a mnemonic for 
the basis in use on a given index, this vertex is accompa¬ 
nied by a normalisation factor and must be manipulated 
in accordance with the rules specified in Refs. 28 and 29. 
We then define the two-index tensor F 3 '^ 7 of Fig. 4(iv) by 
absorbing both the associated numerical coefficients and 
the vertex into c“ [Figs. 4(iii)-(iv)]. Because the charges 
on the vertex absorbed into c ,/ ^ 7 may always be uniquely 
determined from those appearing on the remaining im¬ 
plicit trees, there is a 1:1 mapping between values of a 
and pairs /? 7 , and if we associate 

ol = (o, [l a ) 

P = (&,W>) 

7 = (c,/u) 



FIG. 4. Splitting and fusing of indices: (i) One-index rep¬ 
resentation of a state of five anyons on the disc. The fu¬ 
sion tree basis has been chosen for subsequent convenience, 
(ii) The trunk vertex is made explicit, and (iii) normalisation 
factors are absorbed into a redefinition of the vector, c“ H> c' a . 
(iv) The vertex is then absorbed to yield a two-index repre¬ 
sentation of the state, (v) Application of a three-site operator 
M to state c. (vi) Application of a normalised fusion vertex 
to the greek indices combines two greek indices back into one. 


then for a given correspondence between a and a pair 

Pi: 



Noting that c“ is a state, we have a = I and the above 
expression simplifies by identifying d a = 1 , b = c, and 
thus db = d c . Operator M may now be applied to the 
state as shown in Fig. 4(v), yielding the new state 

d = M%c" Sl . (14) 

The indices appearing on d" may then be recombined. 
Formally, the absorbed vertex is re-expressed as per 


Fig. 4(iii) and is then surmounted by a normalised fusion 
vertex, as per Fig. 4(vi). The loop is evaluated according 
to the rule 



and these factors are absorbed into the redefined vector 
(c"'")“ where they cancel out the vertex normalisation 
factors of [d Q 3 /(G? ai G ? a2 )] 1 ' /4 from Figs. 4(ii) and (vi). In 
practice, the fusion operation taking d" into d"" is al¬ 
ways equivalent to reversing the splitting process. 

The index-splitting process may also be applied to two- 
index objects to yield three-index objects. It is custom- 
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ary not to produce trivalent vertices having three upward 
or three downward legs, as such objects may always be 
rewritten in terms of vertices having at most two legs 
in the same direction and either an additional identity 
charge or a bend. In this paper we extend the same 
convention to the anyonic tensors we employ. We also 
do not create objects with four or more indices, as the 
labelling of the internal vertices becomes ambiguous for 
non-Abelian anyonic models. Such objects are anyway 
unnecessary for the implementation of any tensor net¬ 
work algorithm, as an n-index conventional tensor may 
always be replaced by an anyonic tensor having three or 
less indices, but n leaves, as we shall see in Sec. IIC. 

The index-fusing process may also be applied to three- 
index tensors to obtain two-index tensors, though we note 
that any time a pair of indices are to be combined, either 
both must be upgoing indices, or both must be down¬ 
going indices, to permit their connection with a fusion 
vertex as per Fig. 4(v). Upgoing and downgoing indices 
may be combined, but to do so it is first necessary to con¬ 
vert both into the same type using a bend, as discussed 
in Sec. IIB 4. 

From here on, to reduce the use of primes we will use 
the same label (e.g. c) to represent an object both before 
and after a single-tensor operation. Thus, for example, 
in Fig. 4 the objects c, c', and c" would all be denoted c. 

b. Mixed normalisation: The reader might wonder 
why, in Sec. IIB3a, we have chosen to combine two dif¬ 
ferent normalisation schemes, namely the implicit fusion 
tree normalisation scheme, in which no factors arise from 
the fusion tree diagrams, and the diagrammatic isotopy 
normalisation scheme used in Refs. 28 and 29. The an¬ 
swer is that in doing so, we obtain the best of both worlds: 
Our vectors, matrices, and three-index tensors exhibit 
a natural normalisation where states satisfy c a c' a = 1, 
the identity operator is described by the identity matrix, 
and index contractions may be performed without having 
to evaluate any diagrammatic factors provided the im¬ 
plicit trees match up. However, by switching to the dia¬ 
grammatic isotopy basis on all vertices absorbed into our 
tensors, if we now connect these tensors into a network 
which permits elimination of the (pale grey) implicit fu¬ 
sion trees, then after eliminating these trees we may ma¬ 
nipulate this network using diagrammatic isotopy. The 
same is also true for any subnetwork normalised entirely 
within the diagrammatic isotopy convention. 

Another advantage to using the diagrammatic normal¬ 
isation convention for vertices absorbed into our tensors 
is that we may use vertical bends to convert upper greek 
indices into lower indices, and vice versa. We may in¬ 
clude vertical bends in our tensor networks, and we may 
introduce cup/cap pairs on greek indices at will, which 
is extremely useful when contracting networks. The im¬ 
plementation of vertical bends is described in Sec. IIB 4, 
and their use in tensor networks is discussed in Sec. IIC. 

For consistency, operations should only ever be per¬ 
formed on regions of a network which all adhere to the 
same normalisation convention. Thus, for example, both 


vertices acted on by an F’-move must be normalised ac¬ 
cording to the implicit tree convention, or both must be 
normalised according to the diagrammatic isotopy con¬ 
vention. Also note that the implicit vertex normalisation 
scheme only has meaning in the context of a greek in¬ 
dex, where the labellings of the implicit fusion tree are 
enumerated by that index. Consequently, the only oper¬ 
ation where we can change the normalisation of a vertex 
is during the fusing or splitting of a greek index. 

c. Using fusion and splitting tensors: In Sec. IIB 3 a 
we found it convenient to describe splitting and fusing in 
terms of operations performed on a tensor such as c“ or 
in which context it may be considered a single-tensor 
operation. However, this process may also be understood 
as contraction with a three-index fusing or splitting ten¬ 
sor [with the latter also subsuming factors arising from 
the loop diagram in Fig. 4(v)]. Defining 



we may write 

C “ p = (iV S p lit ) ap cP c“ = (iV fuse ) “ e e-r ( 18 ) 

where N£ b specifies the anyonic fusion rules 

0x6 — N %b c ( 19 ) 

C 

and 8'- jh is 1 if N° b is non-zero and 0 otherwise. Thus N{ use 
and iVgpiit are analogous to Tf use and T sp ii t in Ref. 33, 
and Aspiit corresponds to N in Refs. 19 and 26. 

We note that our definition of Af use differs from that 
in Refs. 19 and 26, where Af use = N sp lit = N, as we are 
employing a tensor-by-tensor approach in the present pa¬ 
per, and thus the factors ( dbd c /d a Y' 2 arising from the 
loop in Fig. 4(v) are absorbed into state c at the time of 
index fusion. We therefore find it convenient to incorpo¬ 
rate these factors into Nf use . In contrast, in Ref. 26 the 
numerical factors arising from loop contraction are evalu¬ 
ated for the tensor network as a whole, rather than on an 
operation-by-operation basis, and thus for fusing as well 
as splitting, only the normalisation factors [d a /{dbd ( f)Y' A 
are incorporated into N. (In Ref. 19 it is noted that the 
factors associated with the diagrammatic isotopy conven¬ 
tion must be accounted for separately, though no explicit 
recipe was given for doing so. Here we give a practical 
recipe for doing so without having to separately evaluate 
the numeric factors associated with labellings of a ten¬ 
sor network, though this comes at the cost of losing the 
equivalence of 7Vf use and A sp ii t .) 

d. Block structure and fusion priority: Although we 
will not be using the fusing and splitting tensor formal¬ 
ism in this paper, there is one important thing which we 
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can learn from them. First, consider that each greek in¬ 
dex a represents a pair of indices (a,/x a ) corresponding 
to charge and degeneracy respectively. As mentioned in 
Sec. II A, when taken in conjunction with the anyonic 
fusion rules, these indices impose a natural block struc¬ 
ture on anyonic tensors as the entries of the tensors may 
be non-zero only where the index charges (supplemented 
by trivial charges where there are less than three indices 
on a tensor) are consistent with a valid fusion tree ver¬ 
tex. For example, given an operator Mp and adhering to 
the convention that all orientation arrows point upwards, 
the vertex and the block structure associated with this 
operator are shown in Table I. When fusing indices, it is 
desirable to explicitly preserve this block structure. Thus 
a vector c ^ 7 will exhibit a block structure with regards 
to the associated charges b and c, having non-zero entries 
only where b = c, but the single-index representation c“ 
should also exhibit block structure with regards to the 
charge a. This is best illustrated with a specific example, 
so let us posit the Z 2 fusion rules with charge labels 0 
and 1 satisfying 


and Hb and /x c taking only a value of 1 . An example map¬ 
ping between pairs /3y and values of a which preserves 
the block structure is given in Table II. This particular 
mapping was constructed as follows: 

1. Let a take the value 0. 

2. Select the first fusion rule in Eq. (20) consistent 
with an output charge of a. 

3. Let b and c be the left and right input charges to 
this rule respectively. 

4. Select the first value of /3 consistent with charge b. 

5. Select the first value of 7 consistent with charge c. 

6 . Associate the pair /3y with the first unassigned 
value of a. 

7. Select the next value of 7 consistent with charge 
c and return to step 6 . If no further compatible 
values of 7 are available, advance to step 8 . 


0 x0^0 
0 x 1 —x 1 
1 x0—>1 
1 x 1 —x 0, 


( 20 ) 


TABLE I. Entries in the matrix representation of an operator 
Mp may be non-zero only where the corresponding charge 
labels yield a valid labelling of a fusion tree vertex. Adhering 
to the convention that all orientation arrows point upwards, 
and inserting a third, trivial charge as M has only two indices, 
we obtain the vertex and block structure associated with M. 
For this example it is assumed that there are only two charges 
in the anyon model, 0 and 1 . 

Vertex structure: 

b 

Block structure: Entries which may be non-zero are marked 
with a star: *. 
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8 . Select the next value of /3 consistent with charge 
b and return to step 5. If no further compatible 
values of /3 are available, advance to step 9. 

9. Select the next fusion rule in Eq. (20) consistent 
with an output charge of a, and return to step 4. If 
no further compatible rules are available, advance 
to step 10 . 

10. Increment the value of a. If this results in an invalid 
value of a, stop. Otherwise return to step 2. 


TABLE II. Fusion of two indices /3 and 7 into a single index 
a requires a mapping between pairs /J 7 and values of a. It is 
preferable that this mapping be made in a manner which pre¬ 
serves the block structure of the tensor, such as the example 
given here. 


a 

(a, fl a ) 

0 ,7 

(6, Mb) 

(c, Me) 

1 

(0,1) 

1,1 

( 0 , 1 ) 

( 0 , 1 ) 

2 

( 0 , 2 ) 

1,2 

( 0 , 1 ) 

(0,2) 

3 

(0,3) 

2,1 

(0,2) 

( 0 , 1 ) 

4 

(0,4) 

2,2 

(0,2) 

(0,2) 

5 

(0,5) 

3,3 

(1,1) 

(1,1) 

6 

(0,6) 

3,4 

(1,1) 

(1,2) 

7 

(0,7) 

4,3 

(1,2) 

(1,1) 

8 

(0,8) 

4,4 

(1,2) 

(1,2) 

9 

(1,1) 

1,3 

( 0 , 1 ) 

(1,1) 

10 

(1,2) 

1,4 

( 0 , 1 ) 

(1,2) 

11 

(1,3) 

2,3 

(0,2) 

(1,1) 

12 

(1,4) 

2,4 

(0,2) 

(1,2) 

13 

(1,5) 

3,1 

(1,1) 

( 0 , 1 ) 

14 

(1,6) 

3,2 

(1,1) 

(0,2) 

15 

(1,7) 

4,1 

(1,2) 

( 0 , 1 ) 

16 

(1,8) 

4,2 

(1,2) 

(0,2) 


P 

(6, Mb) 

1 

( 0 , 1 ) 

2 

(0,2) 

3 

(1,1) 

4 

(1,2) 


7 

(c, Me) 

1 

(0,1) 

2 

(0,2) 

3 

(1,1) 

4 

(1,2) 



































In such a construction, we describe 7 as the fast-cycling 
index and /? as the slow-cycling index. 

Now consider the definition of hermitian conjugation 
given in Sec. II A. Consistency requires that if the fusion 
of upgoing indices is performed using (iVf use )^ 7 , then the 
fusion of downgoing indices is performed using the hermi¬ 
tian conjugate, (Nf use )P 7 . Consequently, if the rightmost 
index of a pair is fast-cycling during fusion of upgoing 
indices, then the rightmost index of a pair must also be 
fast-cycling during fusion of downgoing indices, and this 
applies both to charge (e.g. b) and degeneracy (e.g. mf) 
indices. When fusing indices in software such as Mat- 
lab, it is important to recognise that the software can¬ 
not distinguish between upgoing and downgoing indices. 
Instead, there is a numbered ordering assigned to the in¬ 
dices on a tensor, and in Matlab when combining a pair 
of indices the first index is fast-cycling and the second is 
slow-cycling. Thus it is important that if numbers are 
assigned to upper indices from left to right, then lower 
indices must also be numbered from left to right, in order 
to ensure that fusion is performed in a consistent fash¬ 
ion. One may, of course, number indices in any order if 
appropriate fast and slow cycling is manually enforced. 

Finally, having made the block structure of anyonic 
tensors explicit in Table I, it is useful to point out that 
the total charge index a and its associated degeneracy in¬ 
dex y, a may be identified with indices a and a a in Ref. 18 
below Eq. (4). The third index in this reference, m a , is 
omitted for anyon models. This is because Ref. 18 dis¬ 
cusses the block decomposition of tensors having group 
symmetries, and for symmetry groups this index enumer¬ 
ates states within an irrep, while for many anyon models 
no comparable structure exists. 



FIG. 5. Bending converts a lower index into an upper index. 
Implicit fusion trees are not shown, (i) A tensor T is adjacent 
to a bend vertex. This vertex involves the identity charge and 
the dual to the charge on the bending leg. (ii) Consider the 
fusion diagram associated with each charge labelling of (i) in 
turn, (iii) Since the vertex contained within the tensor is nor¬ 
malised according to the diagrammatic isotopy convention, 
the legs corresponding to greek labels respect the rules of dia¬ 
grammatic isotopy. We may therefore deform the diagram as 
shown, and introduce an additional line carrying the trivial 
charge, (iv) We then perform an F-move, and contract the 
resulting loop using Eq. (15). (v) The resulting trivalent ver¬ 
tex yields the structure of the resulting tensor. The numerical 
coefficients of diagram (iv) are absorbed into the entries of the 
tensor. Note that charge a has been mapped into its dual. 


4- Vertical bending 


The next operation we will consider is vertical bend¬ 
ing of a greek index, converting upper indices into lower 
indices and vice versa. A recipe for evaluating vertical 
bends in fusion tree diagrams is given in Ref. 28 [though 
note there is a typo on the lower index of F in Eq. (2.31)], 
using the notation of Frobenius-Schur indicators. For 
simplicity we shall avoid this notation, adhering instead 
to our convention of labelling the orientations of all line 
segments with upward-pointing arrows, though our treat¬ 
ment will be equivalent. 

a. Evaluation of a bend: Since the indices appear¬ 
ing explicitly on a tensor are normalised according to 
the diagrammatic isotopy convention, if a vertical bend 
appears immediately adjacent to a tensor then we may 
evaluate it as shown in Fig. 5. Indeed, this procedure 
and its clockwise counterpart serve to define the tensors 


in Ref. 28 denoted and [Bf b }^ v according to 



where 


( 21 ) 

( 22 ) 


Ka [-^o ] (I,1,1)(I,1,1) 


(23) 


and >r* appears in Eq. (21) due to the use of a left-facing 
Frobenius-Schur indicator when specifying the diagram 
associated with in Ref. 28. For unitary anyon 

models, Eq. (22) simplifies to the form given in Ref. 28 
because for unitary models, ([ F t '^ bc ]~ 1 )^ = Fff bc . 

Note that, as seen in Fig. 5(i), a bend may be un¬ 
derstood as a vertex with a trivial leg, and this vertex 
must be normalised according to the diagrammatic iso- 
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topy convention and hence is drawn in black. For sim¬ 
plicity, when bending we will generally omit drawing the 
vertex leg carrying the trivial charge. 

b. Introducing new bends: One might hope that the 
approach given in Fig. 5 would suffice for the evalua¬ 
tion of tensor network diagrams which contain bends. In 
practice, however, we often need to introduce additional 
bends when performing tensor contractions, as it is often 
necessary to temporarily lower an index then later raise 
it again. For this we need the identities 



which are valid only if the bending vertices are nor¬ 
malised in the diagrammatic isotopy convention. For now 
we concern ourselves primarily with the introduction of 
bend pairs on an existing greek index. However, in prin¬ 
ciple we may introduce additional bend pairs anywhere 
we like in a fusion diagram (or in Sec. IIC, anywhere 
in an anyonic tensor network), provided we respect the 
normalisation requirement. This is because the insertion 
of a bend pair in the middle of an implicit tree has no 
effect on the number of degrees of freedom involved in 
the labelling of that tree. 

As mentioned in Sec. IIB 3 b, we may apply valid op¬ 
erations to any part of a fusion tree diagram provided 
all vertices within that region share a common normal¬ 
isation. Thus, if we insert a bend pair (normalised in 
the diagrammatic isotopy convention) into a portion of 
the fusion tree which is normalised in the implicit tree 
convention, we may continue to perform our full range of 
operations on parts of the tree adjacent to the inserted 
bends, so long as we do not involve the bending vertices 
themselves. For example, if a greek index on a tensor is 
adjacent to a bend, one may continue to combine that in¬ 
dex with others in the normal way as shown in Fig. 6(i). 
In this Figure, two greek indices on tensor A" are fused 
into a single greek index, and the normalisation conven¬ 
tion of the fusion vertex is changed to implicit. Looking 
exclusively at tensor A" on which the fusion operation 
acts, this is no different to any other fusion operation, as 
shown in Fig. 6(ii). The bend, however, even though it 
is now separated from tensor A ", remains normalised in 
the diagrammatic isotopy convention at all times. If mul¬ 
tiple fusions are performed under a bend, we may even 
perform F-moves on the resulting tree [e.g. Fig. 6(iii)], 
because once again, we may restrict our attention only 
to the portion of the network immediately attached to 
the tensor on which we are operating, giving a unidi¬ 
rectional fusion tree which is normalised entirely in the 
implicit vertex normalisation convention. We shall see 
more worked examples of the insertion of bends when we 
discuss tensor networks in Sec. IIC. 

c. Bending and index ordering: If a greek index car¬ 
ries an implicit fusion tree, note that bending this index 
will cause the branch which was leftmost on the implicit 
tree before bending to be rightmost after bending, and 
vice versa. Consequently, if the labellings of this tree are 



FIG. 6. (i) Example of manipulations involving bends per¬ 
formed on an anyonic tensor network, (ii) Restricting our at¬ 
tention just to tensor T in the second step of diagram (i), the 
fusion and splitting operations are seen to be entirely routine. 
The bend does not participate in these operations, and is no 
more than a deferred instruction to lower the corresponding 
leaf at a later time, when it is once again present on a tensor 
as a greek index, (iii) Similarly, we may perform operations 
such as F-moves and braids so long as all relevant vertices are 
in a consistent normalisation convention. 


ordered according to Sec. IIB3d before bending, then 
they will not be ordered according to this scheme after 
bending, differing by a tree-dependent permutation. We 
avoid explicit calculation of this permutation by recog¬ 
nising that it is always possible to use an appropriate 
combination of F-moves, fusing, and splitting to sepa¬ 
rate the leaves off one at a time and bend them individ¬ 
ually. This achieves the same transformation while only 
bending greek indices corresponding to individual leaves, 
which have no associated implicit fusion tree. 


5. Braiding 

The final single-tensor manipulation we require is the 
braiding of two greek indices. As with bending, this fol¬ 
lows directly from the equivalent manipulation on a fu¬ 
sion tree. Given a tensor if we wish to braid index a 
over index /3 then the entries of acquire coefficients 
based on the values of a in a, b in /3 , and c (and /r, if 
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applicable) in 7 : 



Note that comparison with the equivalent expression for 
fusion trees, 


b a 



c 


shows that /r is the degeneracy index associated with the 
instance of 7 appearing on the left-hand side of the equa¬ 
tion, while v is the degeneracy index associated with the 
instance of 7 on the right. 

Braiding may also be performed within an individual 
implicit fusion tree, and doing so may occasionally sim¬ 
plify the evaluation of a given tensor network. It is never 
obligatory, however, and in the interest of keeping our 
set of basic operations as compact as possible, we shall 
not concern ourselves with it here. 

6. 360° rotation 

Although this identity is never necessary , it can some¬ 
times save computational cost to recognise that when 
bends are applied to the greek indices of a tensor which 
are equivalent to performing a full 360° rotation of that 
tensor, this operation is equivalent to the identity . 6 

C. Tensor networks 

1. Construction 

Having introduced a graphical notation for anyonic 
tensors in Sec. II A, we now introduce anyonic tensor net¬ 
works. To construct an anyonic tensor network, the di¬ 
agrammatic representations of two or more anyonic ten¬ 
sors are drawn with some leaves of their fusion trees con¬ 
nected. Vertical bends may be included, normalised ac¬ 
cording to the diagrammatic isotopy convention. Aside 
from tensors and vertical bends, the only other vertices 
appearing should be those appearing on the unidirec¬ 
tional fusion trees associated with the greek indices of 
the anyonic tensors. Note that it is not necessary to con¬ 
nect all leaves associated with a greek index on a tensor, 
and leaves associated with a greek index may be con¬ 
nected to any number of other tensors. An example of 
an anyonic tensor network is shown in Fig. 7. 



FIG. 7. A simple example of an anyonic tensor network, made 
up of three tensors: w , o, and w. For illustrative purposes, 
in this diagram we have grouped tensors with their associated 
implicit fusion trees using dashed boxes. 

2. Normalisation and diagrammatic isotopy 

In Sec. IIB 3 we adopted a normalisation scheme com¬ 
bining two different vertex normalisations: “Implicit” 
vertices are drawn in grey and their contraction yields no 
numerical factors, giving rise to convenient tensor iden¬ 
tities such as c“ct = 1 for a state vector. “Explicit” 
vertices, such as those absorbed into tensors during fus¬ 
ing and splitting operations, are drawn in black and are 
associated with the factors of the diagrammatic isotopy 
convention . 28 Many operations behave equivalently in 
both normalisations, in particular 

• F -moves on unidirectional fusion trees (all fusion 
vertices or all splitting vertices) with both vertices 
in the same normalisation convention, 

• horizontal deformation of world lines, 

• vertical sliding of diagram portions where this does 
not cause overlap or bending (this operation is triv¬ 
ial), and 

• braiding. 

Although we cautioned earlier against applying opera¬ 
tions to mixed fusion trees, we now note that because of 
this equivalence, horizontal deformation, braiding, and 
vertical sliding operations may be applied to any part 
of a network, regardless of whether normalisation con¬ 
ventions are consistent throughout. By combining the 
horizontal deformation of world lines with braiding, we 
may also horizontally push world lines across tensors. 
(Strictly speaking, E-moves can also be applied to mixed 
trees, but the bookkeeping in the vicinity of bends can 
be confusing, so we caution against this.) 

Regarding operations which differ between the two 
conventions: 
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• The elimination of loops is formally different but 
yields equivalent numerical results, because con¬ 
verting vertices from implicit to explicit normal¬ 
isation introduces numerical factors which exactly 
cancel out those arising from elimination of the loop 
in the diagrammatic isotopy convention. 

• The insertion of additional vertices emitting or ab¬ 
sorbing the trivial charge is not supported in the 
implicit tree convention, as the enumerated la¬ 
bellings do not include this trivial charge. 

• The ability to insert, delete, and evaluate (absorb) 
vertical bends also differs significantly between the 
two conventions. 

In particular, in the implicit tree convention, inconsisten¬ 
cies may arise if an attempt is made to introduce bends, 
connect them to the rest of the fusion tree through ver¬ 
tices involving trivial charges, and then manipulate these 
vertices using f-moves. Consequently we only intro¬ 
duce bends in the diagrammatic isotopy convention, only 
rotate tensors in the diagrammatic isotopy convention, 
and—for simplicity’s sake—never perform operations on 
fusion trees which involve more than one normalisation 
convention. We can, however, continue to apply the full 
range of diagrammatic isotopy operations to sw&diagrams 
which are entirely normalised in the diagrammatic iso¬ 
topy convention. 

Now let us consider the operations available to us from 
the perspective of their actions on the greek indices of an 
individual tensor, with F -moves being considered as a 
change of basis on a single greek index. We see that we 
may: 

• fuse and split indices, 

• use F -moves at trunk level (and also deeper in the 
fusion tree, though this is never actually necessary) 
to affect the manner in which combined indices sep¬ 
arate on splitting, 

• use horizontal deformations and braiding to re¬ 
order greek indices, and 

• use bends to raise and lower greek indices, 

all within a consistent normalisation scheme. In fact, 
when applied in close proximity to the tensor like this, all 
of these operations act on vertices normalised according 
to the diagrammatic isotopy scheme with the exception 
of the trunk-level F-move, which acts on implicit trees. 

Applied appropriately, and supplemented by the 
pairwise contraction operation discussed in Secs. IIA 
and IIC 3, this dictionary of fundamental operations is 
sufficient to contract any anyonic tensor network satis¬ 
fying the prescriptions of Sec. IICl, and we need never 
evaluate operations acting on more remote parts of a ten¬ 
sor’s fusion tree. 



FIG. 8. Example contraction procedure for a pair of tensors 
connected by one leaf. Factors of yc a and Xb are written next 
to the tensors into which they are subsequently absorbed. 


3. Pairwise contraction 

Any pair of tensors may be contracted together into a 
single tensor. In the following, we show how contractions 
are perfomed by means of illustrative examples. 

a. First example: Simple connected pair. Consider 
the network given in Fig. 8 where two tensors are con¬ 
nected by their leaves. This example is slightly more 
complex than the contractions considered in Sec. II A, as 
the tensors in Fig. 8 are connected only by leaves, rather 
than by greek indices. However, using the single-tensor 
manipulations of Sec. IIB, these tensors may be put into 
a form where all summed leaves (and only the summed 
leaves) appear on a single greek index on each tensor, 
and there are no more than three unsummed greek in¬ 
dices in total. The contraction may then be implemented 
by means of a summation over that greek index, with the 
resulting object being a valid anyonic tensor itself having 
no more than three greek indices. The resulting object 
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FIG. 9. COLOUR ONLINE: (i) Naieve construction of an 
anyonic superoperator from w and w' of Fig. 7. (ii) Alterna¬ 
tive construction of the superoperator as a convex tensor S. 
The action of S on an operator o is as shown. In this Figure 
and in Fig. 11, where cups, caps, and factors x originate in 
an application of Eq. (24) we have coloured them in sets to 
show which go together, and assigned labels to the charges on 
the bending indices. The bending vertices continue to be nor¬ 
malised according to the diagrammatic isotopy convention. 
Where we choose to absorb the cup into one tensor and the 
cap into another, we may freely choose whether to absorb the 
corresponding factors x along with the cup or the cap. In 
these examples we have chosen to absorb x with the cap. 

may then be reshaped into any desired form. 

b. Second example: Non-convex network. Next, 
consider the example tensor network given earlier in 
Fig. 7. In this network, we have three tensors: w, w', 
and o. Suppose that we know w and w', and wish to 
find a tensor o such that 

wow t oc o. (27) 

We could trial different operators for o and use a Lanczos- 
type algorithm to obtain those which satisfy Eq. (27), 
first contracting o with w and then the resulting object 
with w', or vice versa, and this would pose no greater 
challenge than did Fig. 8. Alternatively, however, we 
might wish to contract w with w ' to obtain a superoper¬ 
ator, write this in matrix form, compute its eigenvectors, 
and then map these eigenvectors back to operators satis¬ 
fying Eq. (27). If we take this approach, though, then we 
encounter a complication. The naieve construction of the 
superoperator, given in Fig. 9(i), does not seem to corre¬ 
spond to any sort of anyonic tensor we have discussed so 
far. This is because our formalism supports only convex 
tensors. A convex tensor is one where 



FIG. 10. Examples of tensors which are not convex tensors. 

1 . a bounding box may be drawn which completely 
encloses the fusion tree of the tensor without inter¬ 
secting it, 

2 . this bounding box touches the free ends of all leaves 
of the (unidirectional) fusion trees, and 

3. all upgoing leaves are consecutive on this boundary 
and all downgoing leaves are consecutive on this 
boundary. 

Thus tensors having forms such as those shown in Fig. 10 
are excluded. Nevertheless, it is still possible to proceed 
with this calculation. Using the permitted operations 
described in Sec. IIC 2 we may deform this diagram into 
the form of Fig. 9(ii), and an operator o which is an 
eigenoperator of Fig. 9(i) will also be an eigenoperator 
of superoperator S in Fig. 9(ii), where S is defined, and 
acts on o, as shown. 

A more satisfying though slightly more complicated 
approach is shown in Fig. 11 (i). The tensors labelled H 
are identity operators, which may be freely inserted into 
a tensor network at any time, even when (as here) they 
end up separating an inserted bend pair. 

To understand why we have inserted these identity 
operators, let us examine the topmost operator more 
closely. Its lower vertex serves to fuse together two in¬ 
dices during the construction of S'. This then necessarily 
yields a single greek index. The body of the identity op¬ 
erator, 12, is then absorbed onto this greek index to yield 
the upper greek index of S' (though in practice this step 
is trivial as is just the identity matrix). Finally, the 
upper vertex of the identity operator forms part of the 
implicit fusion tree of S'. The insertion of this identity 
operator therefore essentially serves as a prescription as 
to how the upgoing leaves of the network involving w and 
ud are to be fused in order to obtain the upper greek in¬ 
dex of S', and what implicit fusion tree is to be associated 
with this index. 

The other two identity operators behave similarly, with 
one constructing the lower implicit fusion tree of 5", and 
the other the implicit fusion tree of v. As can be seen 
from Fig. ll(ii), the eigenvectors v a corresponding to 
eigenoperators o may then be obtained simply by diag- 
onalising S'p 1 , with operators o ^ recovered from v a as 
per Fig. ll(iii). This example illustrates a useful general 
principle: The anyonic tensor formalism presented here 
and in Ref. 19 supports only convex tensors, but any 
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FIG. 11. COLOUR ONLINE: (i) The superoperator may also 
be written as a matrix S'. If the operator o is rewritten as 
a vector v as shown, then the action of S' on v is simply 
matrix multiplication, S'v, as shown in diagram (ii). Eigen- 
operators o may be recovered from eigenvectors v as shown in 
diagram (iii). 


anyonic tensor calculation can always be re-expressed in 
terms of convex tensors, so there is no loss of generality. 

c. Third example: Disjoint tensors. Finally, we note 
that where a pair of tensors are not connected by any 
leaves, pairwise contraction is achieved by fusing greek 
indices if necessary, so that there are a maximum of two 
greek indices per tensor, inserting an additional greek 
index carrying the trivial charge with no degeneracy onto 
each tensor, upgoing on one tensor and downgoing on the 
other, and then contracting over this index. The indices 
should be inserted such that the resulting tensor will be 
convex. An example is shown in Fig. 12. 




FIG. 12. Example of contraction of a pair of disjoint tensors. 


d. Caveat: It is important to note that where pairs 
of anyonic tensors are contracted by index summation, 
this summation has so far always been over a single greek 
index. When summation is performed over a single greek 
index, the associated fusion tree is free of loops, and so 
no extra numerical factors arise from the diagrammatic 
isotopy convention. In contrast, summing over two greek 
indices at once does create a diagrammatic loop, and thus 
(for example) if A a P is obtained by splitting index 7 on 
A 7 and B a p is obtained by splitting index 7 on then 

A^B,= Y A^B {c , Me) (28) 

C=I, 11 c 

but 

A ap B a p= Y A^A^) B{a ^ )(b ^ )da (29) 

<2, Mo 
b—a, fib 

where the extra factor of d a arises from the loop. For sim¬ 
plicity we therefore advocate always fusing greek indices 
before summation. 
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4- Preferential nature of pairwise contraction 

In Appendix B of Ref. 34 it was shown that any net¬ 
work of ordinary tensors may always be optimally con¬ 
tracted by means of a sequence of pairwise contractions. 
We now argue that the same holds true for anyonic tensor 
networks. 

First, consider that any anyonic tensor admits a de¬ 
composition into labelled fusion tree and degeneracy 
components. In analogy to Eq. (4) of Ref. 18 we may 
write 

¥ “ 0 <g> Y ] (30) 

a,b,c \ c ' 

where V is the space of anyonic tensors, is the de¬ 
generacy space associated with charges a, b , and c, and 
has dimension max (fi a ) x max (/ ib ) x max (/z c ) for a given 
triplet (a,b,c). An anyonic tensor then admits the 
decomposition 

Tf = ® Y (31) 

C 

where for fixed ( a,b,c ) is the degeneracy ten¬ 

sor associated with a given charge sector of T°P, with 
entries enumerated by the values of fi a , fib, and fi c . Ex¬ 
tending this decomposition to the entire tensor network 
then yields a space of labelled fusion diagrams and, for 
each labelling, an associated network of degeneracy ten¬ 
sors. 

A relatively popular approach to the contraction of 
an anyonic tensor network (which we ourselves used in 
Ref. 26) is to iterate over all valid labellings of the fu¬ 
sion tree diagram and, for each labelling, contract the 
associated network of degeneracy tensors using a series 
of pairwise contractions, then multiply by the factor as¬ 
sociated with the labelled fusion tree diagram for the 
entire network. This approach is conceptually simple, 
but is, in general, suboptimal (though for the TEBD al¬ 
gorithm considered in Ref. 26 the overhead is relatively 
small, being at most a factor of the square of the number 
of anyonic charges). Once again, the optimal approach 
is a series of pairwise contractions in which we contract 
together two anyonic tensors at a time, considering only 
the portion of labelled fusion diagram associated with 
this pair and the degeneracy tensors associated with this 
pair. To see how this saving arises, consider the simple 
contraction shown in Fig. 13, and let us nai'evely evaluate 
tensor (AB) without first eliminating the loop. Assuming 
for the moment that no shortcut exists for determining 
valid charge pairs b , c consistent with given values of a 
and d, one must iterate over all possible charges for each 
of b and c, checking consistency at the vertices a, b , c and 
b, c, d. Matters are simplified slightly by applying global 
conservation of charge to show that a = d, but still, if 
the anyon model admits n q charges, for each value of a 



FIG. 13. This contraction provides an example of the im¬ 
portance of eliminating loops at the first available opportu¬ 
nity during tensor contraction (see text), and by implication, 
favouring pairwise contraction of both the degeneracy and the 
fusion tree components of tensors during the evaluation of a 
tensor network. 


one must iterate over ( n q ) 2 labellings of the fusion dia¬ 
gram. If A and B are part of a more extensive network, 
then in general the number of labellings of this network 
which must be checked will scale as 0[(n g ) a ], a > 2. 
Contracting A with B to yield a single tensor (AB) bear¬ 
ing only charge indices a and d thus reduces the number 
of labellings for the entire network by a factor of ( n q ) 2 
to 0[(n 9 ) a-2 ]. This saving significantly reduces compu¬ 
tational cost, and if using a precomputation scheme to 
avoid repeatedly calculating quantities such as the uni¬ 
tary matrices arising from F-moves and braids, 33 then 
this approach also significantly reduces the amount of 
storage space which precomputed data requires. 

This example is a little unrealistic, because the same 
saving can be achieved by identifying iteration over val¬ 
ues of b, fib, and fi c with iteration over fi a , equivalent to 
eliminating the loop diagrammatically. This ceases to be 
universally possible, however, for loops involving the fu¬ 
sion trees of three or more tensors. Consequently, where 
a fusion diagram includes two or more loops involving 
three of more tensors apiece, the use of pairwise contrac¬ 
tions ensures that these loops are eliminated one at a time 
for a saving in the scaling of both computational cost as 
a function of n q , and precomputation storage, again as a 
function of n q . 

For some anyon models, there exist shortcuts which 
may cause the cost savings to be somewhat less pro¬ 
nounced than would be expected from the above exam¬ 
ple. Consider again the nai'eve calculation given above, 
this time specifically for an Abelian anyon model in the 
series. In such a model one may immediately deter¬ 
mine valid values of c given a and b , thus reducing the 
saving to a factor of O (n q ) rather than 0[(n q ) 2 ]. Similar 
cost savings may also be achieved in loops involving three 
or more tensors, but the scaling penalty associated with 
addressing the tensor network as a whole may never be 
entirely eliminated if the network includes two or more 
loops involving three or more tensors. For this reason, a 
pairwise approach to the contraction of tensors is deemed 
preferable. 
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FIG. 14. Graphical representation of tracing two indices on a 
three-index tensor, T' 3 = Tj a . 


It therefore follows that the trace in Fig. 14 may be eval¬ 
uated according to 

t /3 = T 0a O0L = T^’2) {a,fla) da. (32) 

a,fl a 

III. FINITE DMRG ON THE DISC 
A. Matrix Product State Ansatz 


One important question remains: Given the preferen¬ 
tial nature of pairwise tensor contraction, what is the 
optimal sequence of pairwise contractions by which to 
contract a given anyonic tensor network? For normal 
tensors, this problem is known to be NP-hard. 3 ° Normal 
tensors may be considered to be a special case of anyonic 
tensors, with only the vacuum charge, trivial fusion rules, 
and trivial E-moves and braiding, and so it follows that 
this problem is NP-hard for networks of anyonic tensors 
as well. The best known approach for finding optimal 
contraction sequences is a brute-force computer search, 
with the state-of-the-art search for normal tensors being 
the netcon algorithm given in Ref. 34. Manual searches 
relying on human intuition may also be effective. For 
practical purposes we will assume that if netcon identi¬ 
fies a contraction sequence as being optimal for normal 
tensors, then it will also be optimal or near-optimal for 
anyonic tensors. Certainly, it is simple to confirm that 
for DMRG these sequences display the optimal scaling 
of O (D 3 ) in terms of the DMRG refinement parameter 
D (the bond dimension of the MPS), and so they are 
sufficiently close to optimal for practical purposes. 


5. Anyonic tensor trace 

For completeness we now describe one further single¬ 
tensor operation, namely the anyonic tensor trace. To 
evaluate a trace, it is necessary that the traced leaves 
be collected into one upper and one lower greek index, 
and any leaves not participating in the trace be collected 
(using bends if necessary) into a third greek index, which 
may be either upper or lower (see Fig. 14). Note that 
upper leaves are paired with lower leaves, so all leaves 
involved in the trace are connected by a line making a 
pair of consecutive bends. We note that: 

1. After tracing, the untraced third greek index is a 
single index and thus the resulting object has non¬ 
trivial entries only where the charge label in this 
greek index is the identity. 

2. The fusion diagram associated with tracing is there¬ 
fore always a loop, even when a third index is 
present. It is normalised according to the diagram¬ 
matic isotopy convention, and hence is associated 
with a factor of d a . 


1. Matrices vs. tensors 


Two common notations for the Matrix Product State 
(MPS) Ansatz exist. We review these notations briefly 
in the context of a conventional MPS before proceeding 
to the anyonic case. 

Let £ be a lattice of n sites with local dimension d 
and local basis {|y)}, 1 < j < d. In the first nota¬ 
tion, a family of d matrices is associated with each lat¬ 
tice site. At lattice site i we will denote these d matrices 
Af, 1 < j < d. The probability amplitude associated 
with a state \j 1 j 2 • • ■ jn ) on lattice £ is then given by the 
matrix product 




= A [1 U [2] ... A [n_1I A [nI 


n 3 2 


Jn — 1 


(33) 


where the matrix at each site is selected according to 
the basis element \j 1 j 2 • •. j n )■ Strictly speaking, for open 
boundary conditions all matrices A\^ are row vectors and 

all matrices are column vectors so that Cj 1 j 2 ,,,j n is a 
number. 

In the second notation it is recognised that a family of 
matrices, enumerated by an index j, is in fact equivalent 
to a three-index tensor FjF where j is the index selecting 
the matrix in Eq. (33) and a and b are the row and col¬ 
umn indices of this matrix respectively. In this notation, 
Eq. (33) becomes 


c hh-jn 


(rdA (r^Y ...(rA-d 

\ / aj 1 V J bj 2 V 



(34) 

where repeated indices are summed. Again note that rW 
and rN have one fewer index, indicating that for a given 
value of ji, (F^d )aji is a row vector whose entries are 
enumerated by a and for a given value of j n , (fAI)^ is 
a column vector whose entries are enumerated by d. 

As we have developed a formalism for anyonic tensors, 
we find the latter formulation to be better suited to our 
current circumstances. 


2. Graphical representation of the anyonic MPS 

Having established a formalism for anyonic tensor net¬ 
works, we may now write down the anyonic version of the 
Matrix Product State (MPS) Ansatz with open boundary 
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FIG. 15. Anyonic Matrix Product State Ansatz for n sites 
with open boundary conditions. 


conditions. This is shown for open boundary conditions 
in Fig. 15, and differs from the conventional graphical 
representation of the MPS as follows: 

• Because anyonic states are represented by fusion 
trees with all leaves at the top, the open indices of 
the anyonic MPS are upward-pointing rather than 
downward-pointing. 

• There are no horizontal index lines: All lines on a 
fusion tree have an upper and a lower end. 

• All index lines carry orientation arrows (though 
these may be omitted for models where charges are 
self-dual). For our convenience, we choose these all 
to be upward-pointing. 

• On a conventional MPS, the ability of the Ansatz to 
faithfully represent an arbitrary state is constrained 
by the dimensions of the indices connecting the T 
tensors, which are restricted to be no greater than 
some refinement parameter D. For an anyonic MPS 
one must also specify which charges may appear on 
a given index, and divide the total dimension D 
between these charges. 

To perform DMRG, we will also require a graphical 
representation of a Hamiltonian. In the first instance we 
will consider a Hamiltonian which may be written as a 
sum of local terms. For specificity we choose these terms 
to be two-site, having a graphical representation 



We will subsequently extend our treatment to also in¬ 
clude Matrix Product Operator (MPO) Hamiltonians in 
Sec. IIIC. 


B. DMRG algorithm with open boundary 
conditions 

We consider now how the DMRG algorithm may be 
applied to the anyonic MPS of Fig. 15 to compute an 
approximation to the ground state of a Hamiltonian. As 
in conventional DMRG, the anyonic MPS is converged 
towards a representation of the ground state by means of 
a process of variational optimisation. In one step of the 
optimisation process any number of consecutive T tensors 
may be updated simultaneously, though for reasons sub¬ 
sequently made clear we recommend at least two. The 
first optimisation is therefore of Ti and P2, followed by 
r2 and T 3 , then T 3 and T 4 , and so on. Upon reaching 
the end of the chain, the process reverses, first optimising 
r„_i and r n , then r n _2 and r n _i, and so forth. 

To understand the optimisation procedure, let us con¬ 
sider a specific pair, T 4 and Ts, while iterating from left 
to right. We are seeking a choice of T 4 and Ts which min¬ 
imises the value of (i/j\H\ip), represented schematically by 
Fig. 16(i). Deleting T 4 and T 5 from this figure, we first 
seek a tensor v with four leaves which may be inserted in 
their place [Fig. 16(h)], and which maximises the value 
of {%lj\H\ip) subject to the constraint = 1. If we 

reshape v and v' as vectors v a and v' a and the rest of 
the tensor network as a matrix Mg [similar in principle 
to the way that o and S were reshaped as vector and 
matrix in Fig. 11 (i)], and we define A(“ as the matrix ob¬ 
tained when H is replaced by the identity operator, then 
this takes the form of a generalised eigenvalue problem 

Mv = XNv (36) 

where we are looking for the choice of v which minimises 
the value of A for the specified M and N. A little prepa¬ 
ration in getting to this point ensures that the T tensors 
satisfy the identities given in Fig. 16(iii) -we shall see 
how this is enforced when we obtain the updated ver¬ 
sions of T 4 and Ts—and the problem then simplifies to 
an eigenvalue problem 

Mv = Xv (37) 

where we seek the vector v which minimises A. The con¬ 
struction of the product Mv in a manner which preserves 
the structure of v is given in Fig. 17(i), where L 3 is con¬ 
structed as per Fig. 17(ii) and Rq is constructed analo¬ 
gously using tensors from the region of the Ansatz to the 
right of v. 

Obtaining an eigenvector v using the Lanczos method, 
we reshape it into a matrix as shown in Fig. 18(i) and 
then perform a singular value decomposition [Fig. 18(h)] . 
Matrix S is diagonal, containing the positive real singular 
values, and satisfies the normalisation condition 

= 1, (38) 

or equivalently, if the diagonal entries of S in charge sec- 
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FIG. 16. DMRG algorithm part 1: (i) Graphical representation of with the Hamiltonian represented schematically by 

a dotted box. (ii) Tensors T 4 and Ig are replaced by a tensor v, whose structure is shown in the inset, (iii) Careful preparation 
ensures that tensors IT, F 2 , and IT satisfy some useful identities, as do IT,..., F n . 


tor q are denoted s q ,..., s™, 

EE( 4 X = 1 (39) 


where the factor of d q arises from the loop in the diagram, 
as per Secs. IIC 3 d and IIC 5. 

In general, the number of entries in S may now exceed 
the refinement parameter D. Tensor S is therefore trun¬ 
cated, retaining the D non-zero entries which contribute 
the largest quantities to the sum in Eq. (39), and the 
indices on U and V t which connect to S are truncated 
likewise. The truncated matrix S'trunc is then normalised 


according to 


('S'trunc) / 3 —> 


(» 5 trunc )/■ 


(Strunc)Z (Strunc)- 


(40) 


such that Eq. (38) is again satisfied. If D q entries are 
retained in charge sector q , we say that charge sector q is 
of dimension D q on the indices connecting U with S and 
S with yt. 

As we are iterating from left to right, the diagonal 
matrix 5trunc is now absorbed into and the new T 4 
and T 5 are constructed from U and as shown in 
Fig. 18(iii). Construction of T 4 entirely from the unitary 
matrix U ensures that when updating Ts and Tg, the 
identity given in Fig. 16(iii) now extends to T 4 as well. 
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FIG. 17. DMRG algorithm part 2: (i) Construction of the product Mv. The terms have been simplified using the identities of 
Fig. 16(iii). (ii) Construction of operator L 3 in diagram (i). Construction of R§ proceeds analogously on the right, (iii) After 
updating T 4 and Ts, tensor L 4 is constructed in preparation for the update of Ts and T 6 . 
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FIG. 18. (i) Reshaping of v into a matrix, (ii) Singular value 
decomposition of v into U, S, and tV (iii) Construction of 
updated T (3) and F (4) from U and SV'. 


We also note that the index which connected tensor 
U with tensor S now connects the new i ?4 with the new 
IV The dimensions D q of the different charge sectors 
on this index are determined entirely by the relative con¬ 
tributions of the singular values to Eq. (39), and these 
dimensions may be different to those which were exhib¬ 
ited by the old 1?4 and IV This opportunity to optimise 
the dimensions of the different charge sectors as a part 
of the variational update procedure is the reason why we 
recommend updating at least two T-tensors at once. 

In preparation for the next update step, tensor L 4 is 
now computed from L 3 and T 4 as shown in Fig. IT (iii), 
though tensor L 3 is also retained in memory as it will 
be needed again when the next pass, which iterates from 
right to left, updates tensors T 4 and Ts. Similarly, while 
iterating from left to right we employed tensor Rq in 
Fig. 17(i), which was constructed after updating tensors 
T 5 and Tg on the last pass which was performed from 
right to left. 


C. Matrix Product Operator Hamiltonians 


value decompositions, each separating off one site at 
a time to yield structures such as the three-site MPO 
operator shown in Fig. 19(iii). Extrapolation to four or 
more sites is straightforward. 

The advantages of working with MPO operators are 
readily apparent, even for two-site operators. For exam¬ 
ple, when applying matrix M to vector v in Fig. 17 (i), the 
term involving /134 now appears as shown in Fig. 20(i). If 
the dimension of the index labelled du is less than that of 
d x d, where d represents a physical lattice site, then the 
cost of evaluating this diagram is reduced compared to 
Fig. 17 (i). In addition, when computing L 4 the left-hand 
side of Fig. 20(i) (labelled lhs) may be re-used, as shown 
in Fig. 20(ii). This latter advantage becomes more pro¬ 
nounced when working with Hamiltonians whose terms 
span more than two sites, with only a portion of each op¬ 
erator overlapping the sites for which the T-tensors are 
being updated. For example, suppose that the Hamilto¬ 
nian is made up of a sum of three-site MPO terms. One 
may take the region labelled lhs, which describes the left¬ 
most site of a term acting on sites 3 and up, and extend 
it as shown in Fig. 20(iii) to represent the two leftmost 
sites of the same term acting on sites 3 and up during 
the update of tensors Ts and Tg. Not only can we recy¬ 
cle our contraction of the region labelled lhs in Fig. 20(i) 
from one update round to the next, but within a given 
update round we have partially precomputed the contri¬ 
bution of one term in the Hamiltonian to the evaluation 
of Mv, and this precomputed portion remains the same 
from one iteration to the next. 

Of course, sooner or later the final portion of an MPO 
is appended onto an lhs portion, and the object is then 
added into L , which contains all terms in the Hamiltonian 
lying entirely to the left of the P-tensors being updated. 
We may, however, cache the lhs objects for re-use when 
sweeping back from right to left—and needless to say, a 
similar treatment pertains to the portions of Hamiltonian 
terms which extend to the right of the T-tensors being 
updated. 

Depending on the construction of the Hamiltonian, it 
may be preferable to perform an MPO decomposition 
term by term, on collections of terms together, or possi¬ 
bly even on the Hamiltonian as a whole. However, pro¬ 
vided dn <C D the cost of the DMRG algorithm in all 
cases continues to scale as 0(D 3 ). More involved discus¬ 
sions of the role of MPO decompositions in DMRG may 
be found in Refs. 5 and 36. 


We now introduce Matrix Product Operator 
(MPO) decompositions for individual terms in the 
Hamiltonian . 36 Considering first a two-site operator hij 
such as was used in Sec. IIIB, by collecting together 
the left two indices and the right two indices and 
performing a singular value decomposition this operator 
may be rewritten in the form of Fig. 19(i), as shown 
in Fig. 19(ii). This decomposition may be extended to 
larger operators through the use of repeated singular 


D. Periodic boundary conditions 

A form of the DMRG algorithm for non-anyonic sys¬ 
tems with periodic boundary conditions and a computa¬ 
tional cost scaling as O (D 3 ) is given in Ref. 37. Although 
less accurate for a given value of D than finite or infinite 
DMRG with open boundary conditions, this algorithm 
may be useful when seeking to avoid edge effects and we 
therefore describe its adaptation to anyonic systems here. 
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FIG. 19. (i) Matrix Product Operator (MPO) decomposition of a two-site anyonic operator, (ii) Calculation of an example 
MPO decomposition of a two-site anyonic operator, (iii) MPO decomposition of a three-site anyonic operator. 



FIG. 20. (i) Diagram equivalent to the /134 term in Fig. 17(i), with the Hamiltonian term expressed in MPO form. The internal 
bond of the MPO is indicated as being of dimension d,H, while physical sites are of dimension d and MPS bond indices are of 
dimension D. (ii) Calculation of L 4 with MPOs. Note the re-use of the region marked “lhs” in diagram (i). (iii) If working with 
a Hamiltonian having terms involving more than two sites, the “lhs” for the term acting on sites 3 and up during optimisation 
of T 4 and Ts may be extended to obtain the “lhs” for the term acting on sites 3 and up during optimisation of Ts and Te by 
incorporating the next portion of the MPO operator as shown. 


We begin by specifying the construction of the anyonic 
MPS for a ring of anyons on the disc. First, let the 
anyons lie on a ring as shown in Fig. 21 (i), and then, 
for our convenience, let us deform the disc so that all 
of these anyons are brought to the front side of the ring 
[Fig. 21(h)] . An MPS Ansatz for this system of anyons is 


given in Fig. 21 (iii), and it follows from the construction 
specified in Fig. 21 (i)-(ii) that the periodic translation 
operator for this Ansatz appears as shown in Fig. 21(iv). 
The action of the periodic translation operator on the 
MPS Ansatz may be evaluated using the same method 
as was employed in Ref. 12 for the periodic translation 
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T 

View from here 



FIG. 21. (i) Anyons on a circular lattice on the surface of the 
disc, (ii) Adopting the viewpoint shown, deform the disc so 
that all anyons are brought to the front of the ring, (iii) This 
system may be described using a periodic MPS Ansatz, and 
the geometry specified in diagram (ii) reveals that the periodic 
translation operator acts as shown in diagram (iv). (See also 
Ref. 12.) 


operator on the torus. Using this method, we see that 
the action of Fig. 21(iv) on the periodic anyonic MPS 
is to cyclically permute the T-tensors of the MPS, and 
also to multiply each element of the outside tensor [r„ 
in Fig. 21 (iii)] by a phase factor of Rf a where a is the 
corresponding charge on the physical index of the tensor 
for that element. 

Now that we have an anyonic version of the periodic 
MPS Ansatz, extension of the periodic DMRG algorithm 
to anyonic systems is relatively straightforward. At the 
heart of this algorithm is a transfer matrix approach, 
where each term in the Hamiltonian is written in the 
form of an MPO and for each term, each site of the lat¬ 
tice is then associated with a transfer “matrix” having 
one of the forms given in Fig. 22(i). [Strictly speak¬ 
ing, the objects shown in Fig. 22(i) are transfer tensor 
networks, but they may be expressed in matrix form as 
shown in Fig. 22(ii) .] The algorithm exploits the fact 
that when a long chain of identical transfer matrices is 



FIG. 22. (i) Transfer “matrices” in the periodic DMRG al¬ 
gorithm. (Strictly, transfer networks.) (ii) Example showing 
how these anyonic tensor networks can be expressed in matrix 
form through the application of identity operators. 


multiplied together, only a relatively small number of sin¬ 
gular values make a significant contribution to the end 
result. Consequently, one can approximate the transfer 
matrix (which has dimensions ranging between D 2 x D 2 
and D 2 dn x D 2 dn depending on the MPO contribution) 
by a product denoted UdV where UdV has the same di¬ 
mensions as the original transfer matrix, say D 2 x D 2 for 
example, but d is much smaller, perhaps 4x4, and U 
and V are then D 2 x 4 and 4 x D 2 respectively. 38 

Crucial to this calculation is the ability to compute 
the approximate singular value decomposition UdV for 
a cost of O (D 3 ). An algorithm for doing so is provided 
in Ref. 37, and begins by multiplying the transfer matrix 
(which Pippan et al. denote M, but we shall write as 
T) by a D 2 x p (or D 2 dn x p) matrix x where p is the 
number of singular values to be retained. This algorithm 
functions unchanged for transfer networks put in matrix 
form as per Fig. 22, though we note that it is necessary 
to decide how the singular values will be divided among 
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* * 



FIG. 23. Matrix Product State Ansatz for a ring of anyons 
encircling one cycle of the torus (i.e. without twist). For an 
explanation of the * notation, see Ref. 12. 


the different charge sectors when constructing x. 

Finally, we note a couple of practical considerations: In 
the preferred form of the periodic DMRG algorithm de¬ 
scribed at the end of Ref. 37, the periodic MPS is divided 
into thirds and iteration over the T-tensors is performed 
one third at a time. We observe that it is probably easi¬ 
est to use the periodic translation operator to rotate the 
third of the MPS being updated into the centre of the 
Ansatz, though this is not essential. Finally, we observe 
that the expression for Mv [where M is the matrix of the 
eigenvalue problem, as in Eq. (37), and not the transfer 
matrix, as in Ref. 37] involves tensor traces. These are 
evaluated as described in Sec. IIC 5. 


IV. FINITE DMRG ON THE TORUS 

Generalisation from periodic boundary conditions on 
the disc to periodic boundary conditions on the torus 
may be achieved by following the prescriptions for sur¬ 
faces of higher genus given in Ref. 12. We consider here 
the simplest such scenario, where a ring of lattice sites 
encircles the torus around one cycle only. More complex 
two-dimensional scenarios will be considered in forthcom¬ 
ing papers. In this scenario, the Ansatz of Fig. 21 (iii) for 
periodic boundary conditions on the disc is replaced by 
Fig. 23 for anyons on the torus. Once again, one applies 
the periodic DMRG algorithm to a third of the T-tensors 
at a time, and it may be convenient to use the periodic 
translation operator to place those tensors at the centre 
of the Ansatz. The periodic translation operator now 
takes the form 


* 



and its action is again to cyclically permute the lattice 
sites and multiply by a factor Rf a . The key difference 


lies in the definition of the inner product between a bra 
and a ket. For anyons on the disc this is given by Eq. (5), 
with the charge on index a necessarily being I. For a ring 
of anyons on the torus, the corresponding inner product 
is 



with the tensor generalisation giving rise to the normal¬ 
isation condition for a state: 



= 1 , 

where any charges are permitted on the indices provided 
the corresponding vertices within the tensors are consis¬ 
tent with the anyonic fusion rules. [Note that Eqs. (42)- 
(43) have been expressed in a diagrammatic form consis¬ 
tent with Ref. 12, and that to apply this expression to a 
state constructed from T-tensors requires vertical bend¬ 
ing. For reference, if c were constructed by contracting 
a loop of F-tensors then index a = (a,/r a ) would be of 
dimension D , the dimension of the MPS bond, and index 
/? = (b, /ib) would be of dimension d n , the dimension of 
n physical sites.] 


V. EXAMPLE 

We include one simple example to illustrate the capa¬ 
bilities of anyonic DMRG. The Golden Chain' is one of 
the best understood non-Abelian anyonic models, com¬ 
prising a chain of fixed Fibonacci anyons ( t ) subject 
to a nearest-neighbour interaction which favours pair¬ 
wise fusion into either the r channel (ferromagnetic in¬ 
teraction) or the vacuum channel (antiferromagnetic in¬ 
teraction). This model has been studied extensively 
using exact diagonalisation, 7,11,12 Time-Evolving Block 
Decimation, 26 valence bond Monte Carlo, 39 and the 
Multi-Scale Entanglement Renormalisation Ansatz. 19 
DMRG has previously been used to compute the cen¬ 
tral charge, 7,13 and we apply it now to calculation of the 
ground state energy of finite chains. 

Both the ferromagnetic and antiferromagnetic Hamil¬ 
tonians are critical, and thus interest in this model 
has concentrated primarily on the infinite chain. For 
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FIG. 24. COLOUR ONLINE: We define the quantity 
A(E 0 / L) as the difference between the ground state energy 
per site for a given MPS bond dimension D, and the value in 
the large-H limit. In this graph we plot A(Eo/L) as a func¬ 
tion of D for L = 50 ferromagnetic (x) and antiferromagnetic 
(+) Fibonacci anyon chains. We encounter the limit of nu¬ 
merical precision at D = 50 for the antiferromagnetic chain 
and at D = 55 for the ferromagnetic chain. 



FIG. 25. COLOUR ONLINE: Energy per site for ferromag¬ 
netic (x) and antiferromagnetic (+) Fibonacci anyon chains 
as a function of chain length L, with fixed MPS bond dimen¬ 
sion D = 233. The L —» oo limit for the antiferromagnetic 
model is shown as a dashed horizontal line. 


the present paper, however, we are restricted to finite 
systems. In Figs. 24-25 we therefore study the behaviour 
of the ground state energy per site (Eq/L) as a function 
of refinement parameter D and chain length L. In 
Fig. 24 we see that the ground state energy rapidly 
converges as a function of D. The results obtained 
are stable to over fourteen decimal places by D = 55, 
essentially corresponding to the accuracy limit of our 
numerical eigensolver. Meanwhile for fixed D = 233 we 
see that the energy per site continues to show a strong 


dependency on chain length, indicating significant finite 
size effects, especially for the ferromagnetic coupling. 
As we are essentially independent of D in this regime, 
and the energy per site of the antiferromagnetic chain 
smoothly approaches the L —> oo limit for even L , we 
perform a polynomial fit of Eq/L as a function of L _1 to 
extrapolate the energy per site of the antiferromagnetic 
model in the L —> oo limit. 


Model 

Eq/L at L = 50 

Eq/L for L —» oo 

AFM Golden Chain 

-0.752390661 

-0.763932014 


If we compare the antiferromagnetic result with 
the exact value obtained in Ref. 39, 


, 7T 

a = 2 cos — 

5 

Eq d 2 — 4 r°° sech(7rx) 

L 4 d J_ 00 cosh (2x arccos |) — | 

« -0.763932023, 


(44) 


we find agreement to seven significant figures. For com¬ 
parison, the TEBD algorithm on the infinite chain with 
D = 200 is accurate to eight significant figures. 26 We 
therefore see that when energy per site is a smooth func¬ 
tion of the length of chain, applying the DMRG algorithm 
to a series of finite chains and extrapolating to L —> oo is 
capable of yielding ground state energies with accuracy 
comparable to those obtained using infinite TEBD. As 
with conventional DMRG, one major advantage of the 
anyonic DMRG algorithm is the rapidity with which it 
converges: When generating Fig. 25, the ground state 
energies of the simulations were found in every case to 
have converged to a precision of between fourteen and 
fifteen significant figures after only a single iteration of 
the algorithm (one left-to-right sweep and one right-to- 
left sweep), being limited only by the accuracy of the 
numerical eigensolver employed. 

We may also compare this result with the anyonic 
MERA of Ref. 19. This paper analysed the Golden Chain 
using the scale-invariant 3:1 ID MERA with \ = 8, di¬ 
vided between charge sectors as = 3, Xt = 5. The 
computational cost of this simulation is of a similar order 
of magnitude to anyonic DMRG with D = 256, as the 
cost of variationally optimising the 3:1 MERA without 
approximation scales as 0(y 8 ) while the cost of DMRG 
scales as O(-D 3 ). Using the 3:1 anyonic MERA with 
X = [3, 5] we find a ground state energy per site of 


rp MERA 
E o 

^ L—>oo 


-0.7639304 


(45) 


agreeing with the exact result to only five significant fig¬ 
ures and requiring over 31,000 iterations to reach this 
level of precision. Anyonic DMRG therefore provides 
a much faster route to the calculation of ground state 
energies than does the anyonic MERA; the primary ad¬ 
vantage of the MERA calculation is that it gives direct 
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access to the scaling operators and scaling dimensions of 
the critical theory. 

VI. CONCLUSION 

In this paper we have expanded upon the idea of any- 
onic tensors introduced in Ref. 19, introducing a con¬ 
struction for these tensors which avoids evaluating nu¬ 
merical factors associated with whole tensor networks. 
We then used this construction to implement finite 
DMRG algorithms for open and periodic chains of anyons 
on the disc and torus. 

The primary advantages of the DMRG algorithm are 
its speed and accuracy, coupled with the fact that it 
yields an explicit representation of the lowest-energy 
state identified. We then applied the anyonic DMRG 
algorithm to two of the most widely studied non-Abelian 


anyon models, the antiferromagnetic and ferromagnetic 
Golden Chains on the disc, and showed that the finite 
anyonic DMRG algorithm is capable of yielding highly 
accurate results for ground state energies on finite and 
infinite chains (the latter through finite size scaling). The 
accuracy of our results is comparable to that which can 
be obtained using anyonic TEBD and greater than can 
be obtained with an anyonic MERA, and they are ob¬ 
tained at substantially lesser computational cost. Any¬ 
onic DMRG consequently represents a highly effective 
technique for the numerical study of anyonic systems. 
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